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Abstract 

The widely used large-scale diagonalization method using harmonic oscillator basis functions 
(an instance of the Rayleigh-Ritz method, also called a spectral method, configuration-interaction 
method, or "exact diagonalization" method) is systematically analyzed using results for the con- 
vergence of Hermite function series. We apply this theory to a Hamiltonian for a one-dimensional 
model of a quantum dot. The method is shown to converge slowly, and the non-smooth character 
of the interaction potential is identified as the main problem with the chosen basis, while on the 
other hand its important advantages are pointed out. An effective interaction obtained by a simi- 
larity transformation is proposed for improving the convergence of the diagonalization scheme, and 
numerical experiments are performed to demonstrate the improvement. Generalizations to more 
particles and dimensions are discussed. 

PACS numbers: 73.21.La, 71.15.-m, 31.15.Pf 
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I. INTRODUCTION 



Large-scale diagonalization is widely used in many areas of physics, from quantum 
chemistry- to nuclear physics.- It is also routinelyused to obtain spectra of model quantum 
dots, see for example Refs. ^HE 3Il3LlllIl2 l3Il3|l3Jl8|- The method is based on 
a projection of the model Hamiltonian onto a finite-dimensional subspace of the many-body 
Hilbert space in question, hence the method is an instance of the Rayleigh-Ritz method. 19 
Usually, one takes the stance that the many-body Hamiltonian is composed of two parts 
H an d Hi, treating the latter as a perturbation of the former, whose eigenf unctions are 
assumed to be a basis for the Hilbert space. This leads to a matrix diagonalization prob- 
lem, hence the name of the method. As Hi often contains the interaction terms of the 
model, "perturbing" the electronic configuration states of Hq, the method is also called the 
configuration-interaction method. In the limit of an infinite basis, the method is in principle 
exact, and for this reason it is also called "exact diagonalization". Usually, however, this 
method is far from exact, as H\ is rarely a small perturbation (in a sense to be specified in 
Sec. IIIIEp while limited computing resources yield a tight bound on the number of degrees 
of freedom available per particle. 

In this work we provide mathematical convergence criteria for configuration-interaction 
calculations. More specifically, we address this problem in the case where H is a harmonic 
oscillator (or h.o. for short), concentrating on a simple one- dimensional problem. A com- 
mon model for a quantum dot is indeed a perturbed harmonic oscillator, and using h.o. basis 
functions is also a common approach in other fields of many-body physics and partial dif- 
ferential equations settings in general, as it is also known as the Hermite spectral method. 20 
When we in the following refer to the configuration-interaction method, or CI for short, it 
is assumed that a h.o. basis is used. 

Studying a one-dimensional problem may seem unduly restrictive, but will in fact en- 
able us to treat realistic multidimensional problems as well due to the symmetries of the 
harmonic oscillator. Moreover, we choose a worst-case scenario, in which the interaction 
potential decays very slowly. We argue that the nature of the perturbation Hi, i.e., the non- 
smooth character of the Coulomb potential or the trap potential, hampers the convergence 
properties of the method. To circumvent this problem and improve the convergence rate, we 
construct an effective two-body interaction via a similarity transformation. This approach, 
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also using a h.o. basis, is routinely used in nuclear physics,— i22j22 where the interactions are 
of a completely different nature. 

The effective interaction is defined for a smaller space than the original Hilbert space, 
but it reproduces exactly the lowest-lying eigenvalues of the full Hamiltonian. This can 
be accomplished by a technique introduced by Suzuki, Okamoto and collaborators.— i2£i2&2Z 
Approaches based on this philosophy for deriving effective interactions have been used with 
great success in the nuclear many-body problem.— i22j22 For light nuclei it provides bench- 
mark calculations of the same quality as Green's function Monte Carlo methods or other ab 
initio methods. See for example Ref. |28| for an extensive comparison of different methods 
for computing properties of the nucleus 4 He. It was also used in a limited comparative 
study of large-scale diagonalization techniques and stochastic variational methods applied 
to quantum dots. 29 

We demonstrate that this approach to the CI method for quantum dots yields a consid- 
erable improvement to the convergence rate. This has important consequences for studies of 
the time-development of quantum dots with two or more electrons, as reliable calculations 
of the eigenstates are crucial ingredients in studies of coherence. This is of particular im- 
portance in connection with the construction of quantum gates based on quantum dots.— 
Furthermore, the introduction of an effective interaction allows for studies of many-electron 
quantum dots via other many-body methods like resummation schemes such as various cou- 
pled cluster theories as well. As the effective interaction is defined only within the model 
space, systematic and controlled convergence studies of these methods in terms of the size 
of this space is possible. 

The article is organized as follows: In Sec. [Til the model quantum dot Hamiltonian is 
discussed. In Sec. IIHI we discuss the CI method and its numerical properties. Central to 
this section are results concerning the convergence of Hermite function series.— 1 ^ We also 
demonstrate the results with some numerical experiments. 

In Sec. IIVI we discuss the similarity transformation technique of Suzuki and 
collaborators^ 4 ^' 2 ^ 2 ^' and replace the Coulomb term in our CI calculations with this effective 
interaction. We then perform numerical experiments with the new method and discuss the 
results. 

We conclude the article with generalizations to more particles in higher dimensions and 
possible important applications of the new method in Sec. EJ 
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II. ONE-DIMENSIONAL QUANTUM DOTS 



A widely used model for a quantum dot containing N charged fermions is a perturbed 
harmonic oscillator with Hamiltonian 

N 1 

H = E(-2 V ?+2 l|fjl|2 + ^ ) ) 

JV N 

+E E ^(ik-^id, C 1 ) 

j=l k=j+l 

where fj G M 2 , j = 1, . . . ,N are each particle's spatial coordinate, u(r) is a small modification 
of the h.o. potential ||r|| 2 /2, and U(r) is the Coulomb interaction, viz, U(r) = X/r. Modelling 
the quantum dot geometry by a perturbed harmonic oscillator is justified by self-consistent 
calculations,— 1 ^ 1 ^ and is the stance taken by many other authors using the large-scale 
diagonalization technique as well.-^&^^^iM^ii&il 

Electronic structure calculations amount to finding eigenpairs (E,m), e.g., the ground 
state energy and wave function, such that 

Hm = EV, m G H and E G R. 

Here, even though the Hamiltonian only contains spatial coordinates, the eigenfunction 
\& is a function of both the spatial coordinates G M 2 and the spin degrees of freedom 
a k G {-1/2, +1/2}, i.e., 

H = L 2 (R 2N )®C 2 . 

The actual Hilbert space is the space of the antisymmetric functions, i.e., functions m for 
which 

m(xp {1) ,xp {2 ), ■ ■ -,Xp(N)) = sga(P)m(xi,x 2 , . . .,x N ), 

for all permutations P of iV symbols. Here, x k = (^c? 

For simplicity, we concentrate on one-dimensional quantum dots. Even though this is 
not an accurate model for real quantum dots, it offers several conceptual and numerical 
advantages. Firstly, the symmetries of the harmonic oscillator makes the numerical prop- 
erties of the configuration-interaction method of this system very similar to a two or even 
three-dimensional model, as the analysis extends almost directly through tensor products. 
Secondly, we may investigate many-body effects for moderate particle numbers N while still 
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allowing a sufficient number of h.o. basis functions for unambiguously addressing accuracy 
and convergence issues in numerical experiments. 

In this article, we further focus on two-particle quantum dots. Incidentally, for the two- 
particle case one can show that the Hilbert space of anti-symmetric functions is spanned by 
functions on the form 

*(ri,0-i,f 2 ,0-2) = i>(f 1 ,f 2 )x((Ti,(T2), 

where the spin wave function x can be taken as symmetric or antisymmetric with respect 
to particle exchange, leading to an antisymmetric or symmetric spatial wave function if), 
respectively. Inclusion of a magnetic field B poses no additional complications—, but for 
simplicity we presently omit it. Thus, it is sufficient to consider the spatial problem and 
produce properly symmetrized wavefunctions. 

Due to the peculiarities of the bare Coulomb potential in one dimension~£^ we choose a 
screened approximation U(x\ — x 2 ; A, 5) given by 



where A is the strength of the interaction and 5 > is a screening parameter which can be 
interpreted as the width of the wave function orthogonal to the axis of motion. This choice 
is made since it is non-smooth, like the bare Coulomb potential in two and three dimensions. 
The total Hamiltonian then reads 



Observe that for U = 0, i.e., A = 0, the Hamiltonian is separable. The eigenfunctions of 
H (disregarding proper symmetrization due to the Pauli principle) become ^1(^1)^2(^2), 
where ip n (x) are the eigenfunctions of the trap Hamiltonian H t given by 



Similarly, for a vanishing trap modification v(x) = the Hamiltonian is separable in (nor- 
malized) centre-of-mass coordinates given by 




(2) 




(3) 



X 



Xi + x 2 



and x 



x\ - x 2 



V2 



V2 
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Indeed, any orthogonal coordinate change leaves the h.o. Hamiltonian invariant (see Sec. lIIII) . 
and hence 

v((X + x)/V2) +v((X- x)/V2) + U(V2x; A, 5). 

The eigenfunctions become <p n (X)ip m (z), where (j) n {X) are the Hermite functions, i.e., the 
eigenf unctions of the h.o. Hamiltonian (see Sec. IIIip . and where i/j m (x) are the eigenfunctions 
of the interaction Hamiltonian, viz, 

Hi = -\^ + \j + U{V2x-,\,S). (4) 

Odd (even) functions ip m (x) yield antisymmetric (symmetric) wave functions with respect 
to particle interchange. 

III. CONFIGURATION-INTERACTION METHOD 
A. The harmonic oscillator and model spaces 

The configuration-interaction method is an instance of the Rayleigh-Ritz method,— em- 
ploying eigenfunctions of the unperturbed h.o. Hamiltonian as basis for a finite dimensional 
Hilbert space V, called the model space, onto which the Hamiltonian ([1]), or in our sim- 
plified case, the Hamiltonian ([2]), is projected and then diagonalized. As mentioned in the 
Introduction, this method is in principle exact, if the basis is large enough. 

We write the iV-body Hamiltonian (CQ) as 

H = H + Hi, 
with Ho being the h.o. Hamiltonian, viz, 

N N 
i=l 3=1 

and H\ being a perturbation of Hq, viz, 

N N N 

#i = E^') + E E U(\\ri-r*W)- 

3=1 3=1 k=j+l 
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For a simple one-dimensional model of two particles we obtain 

Hq = h(xi) + h(x 2 ), 

where h(x) is the well-known one-dimensional harmonic oscillator Hamiltonian, viz, 

1 d2 1 2 

h{x) = + r ■ 

Clearly, Hq is just a two-dimensional h.o. Hamiltonian, if we disregard symmetrization due 
to the Pauli principle. For the perturbation, we have 

A 



Hi = v(xi) + v(x 2 ) + 



\xi — x 2 \ + 5 

In order to do a more general treatment, let us recall some basic facts about the harmonic 
oscillator. 

If we consider a single particle in D-dimensional space, it is clear that the D-dimensional 
harmonic oscillator Hamiltonian is the sum of one- dimensional h.o. Hamiltonians for each 
Euclidean coordinate, viz, 

h iD) = + 1 -\\x\\> = X>(**)- (5) 

k=l 

We indicate the variables on which the operators depend by parenthesis if there is danger 
of confusion. Moreover, the h.o. Hamiltonian for iV (distinguishable) particles in d dimen- 
sions is simply hS Nd K The D-dimensional h.o. Hamiltonian is manifestly separable, and the 
eigenfunctions are 

D 

fc=i 

with energies 

D ° 
£rt = — + 2_^n k , 

k=l 

where n denotes the multi-index of quantum numbers n^. The one- dimensional h.o. eigen- 
functions are given by 

<j) n {x) = {2 n n\^ 2 Y Xl2 E n {x)e- x2 l\ 



where H n (x) are the usual Hermite polynomials. These functions are the Hermite functions 
and are treated in further detail in Sec. IIII CI 
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As for the discretization of the Hilbert space, we employ a so-called energy-cut model space 
V, defined by the span of all h.o. eigenfunctions with energies up to a given e = N max + D/2, 
viz, 

V := sp |0<5> fc < iV max }, 

k 

where we bear in mind that the D = Nd dimensions are distributed among the N particles. 
For the one- dimensional model with only one particle, the model space reduces to 



V l = sp 



< n < iV max }. 



(6) 



Thus, one particle is associated with one integer quantum number n, denoting the "shell 
number where the particle resides" , in typical terms. For two particles, we get 

V 2 = sp {0 ni (xi)0 n2 (x 2 ) | < m + n 2 < iV max }. 

We illustrate this space in Fig. [TJ 




FIG. 1: Two-body model space defined by a cut in energy. The two-body state has quantum 
numbers n\ and n 2 , the sum of which does not exceed iV max . 



Proper symmetrization must also be applied. However, the Hamiltonian ([T]) commutes 
with particle permutations, meaning that the eigenfunctions will be symmetric or antisym- 
metric, assuming that the eigenvalues are distinct. In the case of degeneracy, we may simply 
produce (anti) symmetric eigenfunctions by taking linear combinations. 

We mention that other model spaces can also be used; most common is perhaps the 
direct product model space, defined by N direct products of V\ rather than a cut in energy 
as above. 
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B. Coordinate changes and the h.o. 



It is obvious that any orthogonal coordinate change y = Sx where S T S = 1 commutes 
with h^ D \ In particular, energy is conserved under the coordinate change. Therefore, the 
eigenfunctions of the transformed Hamiltonian will be a linear combination of the original 
eigenfunctions of the same energy, viz, 

where the sum is over all n' such that t^> = eg. Here, T performs the coordinate change, 
viz, 

f$ n (x) = (7) 

where T is unitary. Also note that energy conservation implies that V is invariant with 
respect to the coordinate change, implying that the CI method is equivalent in the two 
coordinate systems. 

An important example is the centre-of-mass transformation introduced in Sec. [Til This 
transformation is essential when we want to compute the Hamiltonian matrix since the 
interaction is given in terms of these coordinates. 

Observe that in the case when the Hamiltonian is in fact separated by such a coordinate 
change, the formulation of the exact problem using h.o. basis is equivalent to two one-particle 
problems using h.o. basis in the new coordinates. 

C. Approximation properties of the Hermite functions 

In order to understand the accuracy of the CI method, we need to study the approx- 
imation properties of the Hermite functions. Note that all the Hermite functions <p n (x) 
spanning L2OR) are smooth. Indeed, they are holomorphic in the entire complex plane. Any 
finite linear combination of these will yield another holomorphic function, so any non-smooth 
function will be badly approximated. This simple fact is sadly neglected in the configuration- 
interaction literature, and we choose to stress it here: Even though the Hermite functions 
are simple to compute and deal with, arising in a natural way from the consideration of a 
perturbation of the h.o. and obeying a wealth of beautiful relations, they are not very well 
suited for computation of functions whose smoothness is less than infinitely differentiable, or 
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whose decay behaviour for large |x| is algebraic, i.e., f(x) = o(\x\P) for some (5 < 0. Due to 
the direct product nature of the iV-body basis functions, it is clear that these considerations 
are general, and not restricted to the one-dimensional one-particle situation. 

Consider an expansion i/j(x) = Y2^=o c n4>n{x) in Hermite functions of an arbitrary ip 6 
L 2 (1R). The coefficients are given by 

/■oo 

c n = ((f) n ,ip) = / ^(x)H n (x)e~ x2/2 dx. 



Here, H n (x) = {2 n n\^y x l 2 H n (x) are the normalized Hermite polynomials. If ip(x) is well 
approximated by the basis, the coefficients c n will decay quickly with increasing n. The least 
rate of convergence is a direct consequence of 

oo 

IHI 2 = $>„| 2 < oo, 

ra=0 

hence we must have |c n | = o(n^ 1 ^ 2 ). (This is not a sufficient condition, however.) With 
further restrictions on the behaviour of i[>(x), the decay will be faster. This is analogous 
to the faster decay of Fourier coefficients for smoother functions— although for Hermite 
functions, smoothness is not the only parameter as we consider an infinite domain. In this 
case, another equally important feature is the decay of ip(x) as \x\ grows, which is intuitively 
clear given that all the Hermite functions decay as exp(— x 2 /2). 



Let us prove this assertion. We give here a simple argument due to Boyd (Ref. 1311 ). but 
we strengthen the result somewhat. 

To this end, assume that ip(x) has k square integrable derivatives (in the weak sense) and 
that x m ip(x) is square integrable for m = 0, 1, . . . , k. Note that this is a sufficient condition 
for 

I 

a)il){x) = —=(xil){x) — vjj'(x)), 
v 2 

and (a^) 2 i/j(x) up to {a^) k i}j{x) to be square integrable as well. Here, a) and its Hermitian 
conjugate a are the well-known ladder operators for the harmonic oscillator. 40 
Using integration by parts, the formula for c n becomes 

i){x)H n {x)e- x2/2 dx 

) 

poo 

{n + 1)- 1/2 / [a^(x)}H n+1 (x)e- x2/2 dx, 
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or 

c„ = -{n + iy 1/2 d n+1 , 

where d n are the Hermite expansion coefficients of a*i[)(x) G L2. Since ^2 \dn\ 2 < 00 by 
assumption, we obtain 

00 

n|c n | 2 < OO, 

n=0 

implying 

c n = o(n _1 ). 

Repeating this argument k times, we obtain the estimate 

Cn = ( n ~(k+D/2y 

It is clear that if ip(x) is infinitely differentiable, and if in addition ip(x) decays faster 
than any power of x, such as for example exponentially decaying functions, or functions 
behaving like exp(— ax 2 ), c n will decay faster than any power of 1/n, so-called "infinite- 



order convercence," or "spectral convergence." Indeed, Hille (Ref. |32| ) gives results for the 
decay of the Hermite coefficients for a wide class of functions. The most important for our 
application being the following: If ip(x) decays as exp(— ax 2 ), with a > 0, and if r > is 
the distance from the real axis to the nearest pole of i/j(x) (when considered as a complex 
function), then 

| Cn | = 0(n- 1 / 4 e- T ^+ T ), (8) 

a very rapid decay for even moderate r. 

An extremely useful property^ of the Hermite functions is the fact that they are uniformly 
bounded, viz, 

\<f> n (x)\ < 0.816, Vx,n. 

As a consequence, the pointwise error in a truncated series is almost everywhere bounded 
by 

iVmax OO 

\^{X) - Cn0n(x)| < 0.816 l C "l" 

n=0 ra =A r max +l 

Thus, estimating the error in the expansion amounts to estimating the sum of the neglected 
coefficients. If |c n | = o(n a ), 

ax 



IV'O) - c nM%)\ = o(N^), a.e. 



n=0 
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For the error in the mean, 

N 

\Mx)-Y^c n M*)\\=0{N£Z'>), (9) 

n=0 

as is seen by approximating Y^n=N +1 l c «| 2 by an integral. 

In the above, "almost everywhere", or "a.e." for short, refers to the fact that we do not 
distinguish between square integrable functions that differ on a point set of Lebesgue measure 
zero. Moreover, there is a subtle distinction between the notations 0(g(n)) and o(g(n)). For 
a given function /, f{n) = o(g(n)) if lim^oo \f(n)/g(n)\ = 0, while f(n) = 0(g(n)) if we 
have lim^oo \f(n)/g(n)\ < oo; a slightly weaker statement. 



D. Application to the interaction potential 

Let us apply the above results to the eigenproblem for a perturbed one-dimensional 
harmonic oscillator, i.e., 

^"{x) = [x 2 + 2f(x) - 2£]<0(z), (10) 

which is also applicable when the two-particle Hamiltonian (j2J) is separable, i.e., when U — 
or v — 0. 

It is now clear that under the assumption that f(x) is k times different iable (in the 
weak sense), and that f(x) = o(|x| 2 ) as \x\ — > oo, the eigenfunctions will be k + 2 times 
(weakly) differentiable and decay as exp(— x 2 /2) for large \x\. Hence, the Hermite expansion 
coefficients of ip(x) will decay as o(n a ), a = — (k + 3)/2. 

If we further assume that f(x) is analytic in a strip of width r > around the real axis, 
the same will be true for ip(x), such that we can use Eq. (jSJ) to estimate the coefficients. 

A word of caution is however at its place. Although we have argued that if a given 
function can be differentiated k times (in the weak sense) then the coefficients decay as 
o(n a ), a = —(k + l)/2, it may happen that this decay "kicks in" too late to be observable 
in practical circumstances. 

Consider for example the following function: 

e -xV2 

\x\ + 

which has exactly one (almost everywhere continuous) derivative and decays as exp(— x 2 /2). 
However, the derivative is seen to have a jump discontinuity of magnitude 2/5 2 at x = 0. 
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Coefficients of f(x) for 5 e [0.01, 2] Decay rate a for <5 6 [0.01, 2] 




lnio(n) 



FIG. 2: (Left) Coefficients \c n \ of the function exp(-x 2 /2)/(M + 5) for 5 G [0.01,2], n = 

0. 2. ... , 5000. (Right) Estimated decay rate a of \c n \, i.e., slope of the graphs in the left panel. 

From the theory, we expect o(n~ l ) decay of the coefficients, but for small 5 the first derivative 
is badly approximated, so we expect to observe only o(n -1 / 2 ) decay for moderate n, due to 
the fact that the rate of decay of the coefficients of g(x) are explicitely given in terms of the 
coefficients of a)g{x). 

In Fig. [2] the decay rates at different n and for various 5 are displayed. The decay rate 
a is computed by estimating the slope of the graph of In |c n | versus Inn, a technique used 
thoughout this article. Indeed, for small 5 we observe only a ~ —1/2 convergence in practical 
settings, where n is moderate, while larger 5 gives a ~ — 1 even for small n. 

The above function was chosen due to its relation to the interaction Hamiltonian (jlj). 
Indeed, its coefficients are given by 

Cn = (4>n,g) = (4>n,U(x; 1,5)0 O ), 

1. e., the proportional to the first row of the interaction matrix. Moreover, due to Eq. ( fTUi) . the 
ground state ip of the interaction Hamiltonian has a second derivative with similar behaviour 
near x = as g(x). Thus, we expect to observe a ~ —3/2, rather than a ~ —2, for the 
available range of n in the large-scale diagonalization experiments. 

We remark here, that it is quite common to model quantum dot systems using non- 
smooth potentials^ 1 v (f), and even to use the CI method with h.o. basis functions on these 
models^^ 
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FIG. 3: Left: Symmetric double-well potential created with the Gaussian perturbation 
Aexp[— C(x — fi) 2 ] with A = 4, fj, = and C = 2. Right: Asymmetric double-well potential 
created with the Gaussian perturbation with A = 4, [i = 0.75 and C = 2, and single- well potential 
using C = 0.25. 

E. Numerical experiments 

We wish to apply the above analysis by considering the model Hamiltonian ([2]). We first 
consider the case where v(x) — or U(x) = 0, respectively, which reduces the two-particle 
problem to one-dimensional problems through separation of variables, i.e., the diagonaliza- 
tion of the trap Hamiltonian H t and the interaction Hamiltonian Hi in Eqs. ([3]) and (j4j). 
Then we turn to the complete non-separable problem. 

For simplicity we consider the trap x 2 /2 + v(x) with 

v(x) = Ae- C(x ~^ 2 , A, C > 0, (jt e R, 

which gives rise to a double-well potential or a single-well potential, depending on the 
parameters, as depicted in Fig. [3J The perturbation is everywhere analytic and rapidly 
decaying. This indicates that the corresponding configuration-interaction energies and wave 
functions also should converge rapidly. In the below numerical experiments, we use A = 4, 
C = 2 and \x = 0.75, creating the asymmetric double well in Fig. [3j 

For the interaction Hamiltonian H\ and its potential x 2 /2 + U(y/2x; X,5) we arbitrarily 
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choose A = 1 and 5 = 0.01, giving a moderate jump discontinuity in the derivative. 

As these problems are both one-dimensional, the model space reduces to V\ as given in 
Eq. ($5$). Each problem then amounts to diagonalizing a matrix H with elements 

1 f°° 

H n , m = ((f) n , H tji (f> m ) = (n + -)S n>m + / 4> n (x)f(x)(f> m (x) dx, < n,m < iV max , 

^ J -oo 

with f(x) = v(x) or f(x) = U(y/2x; 1,0.01). We compute the matrix to desired precision 
using Gauss-Legendre quadrature. In order to obtain reference eigenfunctions and eigenval- 
ues we use a constant reference potential method^ implemented in the Matslise package^ 
for Matlab. This yields results accurate to about 14 significant digits. 

In Fig. H] (left) the magnitude of the coefficients of the exact ground states alongside the 
ground state energy error and wave function error (right) are graphed for each Hamiltonian, 
using successively larger iV max . The coefficients of the exact ground states decay according to 
expectations, as we clearly have spectral convergence for the H t ground state, and o(n~ 157 ) 
convergence for the H\ ground state. 

These aspects are clearly reflected in the CI calculations. Both the H t ground state energy 
and wave function converge spectrally with increasing iV max , while for Hi we clearly have 
algebraic convergence. Note that for H t , A^ max ~ 40 yields a ground state energy accurate 
to ~ 10~ 10 , and that such precision would require iV max ~ 10 12 for Hi, which converges only 
algebraically. 

Intuitively, these results are easy to understand: For the trap Hamiltonian a modest value 
of -/V max produces almost exact results, since the exact ground state has extremely small 
components outside the model space. This is not possible for the interaction Hamiltonian, 
whose exact ground state is poorly approximated in the model space alone. 

If we consider the complete Hamiltonian (T5]), we now expect the error to be dominated by 
the low-order convergence of the interaction Hamiltonian eigenproblem. Fig. H] also shows 
the error in the ground state energy for the corresponding two-particle calculation, and the 
error is indeed seen to behave identically to the H\ ground state energy error. (That the 
energy error curve is almost on top of the error in the wave function for H[ is merely a 
coincidence.) 

It is clear that the non-smooth character of the potential U destroys the convergence of 
the method. The eigenfunctions will be non-smooth, while the basis functions are all very 
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Error in ground state energy and wave function 
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FIG. 4: Left: Coefficients of the exact ground states of the Hamiltonians -fft.i- For H{ only 
even-numbered coefficients are nonzero and thus displayed. The almost straight line indicates 
approximately o(n~ 157 ) decay of the coefficients around n = 600 and o(n -1 ' 73 ) around n = 5000. 
Compare with Fig. For the H t ground state we clearly have spectral convergence. Right: 
The error in the ground state energies and wave functions when using the CI method. For Hi 
we have o(n -L24 ) decay for the energy error, and o(n -1 ' 20 ) decay for the wave function error, 
both evaluated at n = 600. For H t we clearly have spectral convergence. A full two-particle CI 
calculation is superimposed, showing that the error in the interaction part of the Hamiltonian ([5]) 
completely dominates. Here, the error in the energy is o(iV~;[ x 02 ) at n = 70, while for H\ alone, we 
have o{N~l^). 



smooth. Of course, a non-smooth potential v(x) would destroy the convergence as well. 

In this sense, we speak of a "small perturbation Hi' if the eigenvalues and eigenfunctions 
of the total Hamiltonian converge spectrally. Otherwise, the perturbation is so strong that 
the very smoothness property of the eigenfunctions vanish. In our case, even for arbitrary 
small interaction strengths A, the eigenfunctions are non-smooth, so that the interaction is 
never small in the sense denned here. On the other hand, the trap modification v (x) repre- 
sents a small perturbation of the harmonic oscillator if it is smooth and rapidly decaying. 
This points to the basic deficiency of the choice of h.o. basis functions: They do not capture 
the properties of the eigenfunctions. 

We could overcome this problem by choosing a different set of basis functions for the 
Hilbert space, and thus a different model space V altogether. However, the symmetries 
of the h.o. lets us treat the interaction potential with ease by explicitly performing the 
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centre-of-mass transformation, a significant advantage in many-body calculations. In our 
one-dimensional case, we could replace U{x\ — X2) by a smooth potential; after all U is just 
an approximation somewhat randomly chosen. We would then obtain much better results 
with the CI method. However, we are not willing to trade the bare Coulomb interaction in 
two (or even three) dimensions for an approximation. After all we know that the singular 
and long-range nature of the interaction is essential. 

We therefore propose to use effective interaction theory known from many-body physics 
to improve the accuray of CI calculations for quantum dots. This replaces the matrix in 
the h.o. basis of the interaction term with an approximation, giving exact eigenvalues in the 
case of no trap perturbation v(x), regardless of the energy cut parameter iV max . We cannot 
hope to gain spectral convergence; the eigenfunctions are still non-smooth. However, we can 
increase the algebraic convergence considerably by modifying the interaction matrix for the 
given model space. This is explained in detail in the next section. 

IV. EFFECTIVE HAMILTON! AN THEORY 
A. Similarity transformation approach 

The theories of effective interactions have been, and still are, vital ingredients in many- 
body physics, from quantum chemistry to nuclear physics: l n fields like nuclear 
physics, due to the complicated nature of the nuclear interactions, no exact spatial potential 
exists for the interactions between nucleons. Computation of the matrix elements of the 
many-body Hamiltonian then amounts to computing, for example, successively complicated 
Feynman diagrams,— 1 ^ motivating realistic yet tractable approximations such as effective 
two-body interactions. These effective interactions are in turn used as starting points for 
diagonalization calculations in selected model spaces.-^^ 1 ^ Alternatively, they can be used 
as starting point for the resummation of selected many-body correlations such as in coupled- 
cluster theories- In our case, it is the so-called curse of dimensionality that makes a direct 
approach unfeasible: The number of h.o. states needed to generate accurate energies and 
wave functions grows exponentially with the number of particles in the system. Indeed, the 
dimension of V grows as N^ x / (Nd)\ 

For the derivation of the effective interaction, we consider the Hamiltonian (j2J) in centre- 
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of-mass coordinates, i.e., 

H = h(X) + h(x) + v((X + x)/V2) + 
v((X-x)/V2) +U(V2x;X,5). 

For v(x) 7^ 0, the Hamiltonian is clearly not separable. The idea is then to treat v(xj) as 
perturbations of a system separable in centre-of-mass coordinates; after all the trap potential 
is assumed to be smooth. This new unperturbed Hamiltonian reads 

H' = h(X) + h(x) + V, 

where V = U(y/2x; A, 5), or any other interaction in a more general setting. We wish to 
replace the CI matrix of H' with a different matrix H' eS , having the exact eigenvalues of H', 
but necessarily only approximate eigenvectors. 

The effective Hamiltonian H' c{{ can be viewed as an operator acting in the model space 
while embodying information about the original interaction in the complete space 7i. We 
know that this otherwise neglected part of Hilbert space is very important if V is not small. 
Thus, the first ingredient is the splitting of the Hilbert space into the model space V = PH. 
and the excluded space Q = QH = (1 — P)H. Here, P is the orthogonal projector onto the 
model space. 

In the following, we let iV be the dimension of the model space V. There should be no 
danger of confusion with the number of particles iV = 2, as this is now fixed. Moreover, we 
let {<&n\n=i be an orthonormal basis for V, and {$n}^ = N+i be an orthonormal basis for Q. 

The second ingredient is a decoupling operator u. It is an operator defined by the prop- 
erties 

Puj = u)Q = 0, 

which essentially means that w is a mapping from the model space to the excluded space. 
Indeed, 

u = (P + Q)u(P + Q) = PuoP + PuoQ + QuoP + QuoQ 
= QluP, 

which shows that the kernel of u> includes Q, while the range of uj excludes V, i.e., that uj 
acts only on states in V and yields only states in Q. 
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The effective Hamiltonian H e g = P[h(x) + h(X)]P + V e g, where V e e is the effective 
interaction, is given by the similarity transformational 

H cS = Pe~ z He z P, (11) 

where z = artanh(c<j — u;*). The key point is that e z is a unitary operator with (e 2 ) -1 = e~ z , 
so that the N eigenvalues of H' eS are actually eigenvalues of H' . 

In order to generate a well-defined effective Hamiltonian, we must define u = QuoP 
properly. The approach of Suzuki and collaborators^ 1 ^ 1 ^ 1 ^ is simple: Select an orthonormal 
set of vectors {Xn}n=v These can be some eigenvectors of H' we wish to include. Assume 
that {PXn\n=i i s a basis for the model space, i.e., that for any n < N we can write 

N 
m=l 

for some constants a n , m . We then define ui by 

UjPXn-=QXn, n=l,...,N. 

Observe that u defined in this way is an operator that reconstructs the excluded space 
components of Xn given its model space components, thereby indeed embodying information 
about the Hamiltonian acting on the excluded space. 

Using the decoupling properties of u we quickly calculate 

N 

UJ® n = QujP^ n = QuJ^ a n,mXm, U = 1, . . . , N 

m=l 

and hence for any n' > N we have 

N 
m=l 

yielding all the non-zero matrix elements of u>. 

As for the vectors Xn, we do not know a priori the exact eigenfunctions of H', of course. 
Hence, we cannot find H' eS exactly. The usual way to find the eigenvalues is to solve a 
much larger problem with N' > N and then assume that these eigenvalues are "exact". 
The reason why this is possible at all is that our Hamiltonian H' is separable, and therefore 
easier to solve. However, we have seen that this is a bad method: Indeed, one needs a 
matrix dimension of about 10 10 to obtain about 10 significant digits. Therefore we instead 
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reuse the aforementioned constant reference potential method to obtain eigenfunctions and 
eigenvectors accurate to machine precision. 

Which eigenvectors of H' do we wish to include? Intuitively, the first choice would be the 
lowest N eigenvectors. However, simply ordering the eigenvalues "by value" is not what we 
want here. Observe that H' is block diagonal, and that the model space contains iV max + 1 
blocks of sizes 1 through iVmax + 1. If we look at the exact eigenvalues, we know that they 
have the structure 

E n , m = {n + l/2) + e m , 

where n is the block number and e m are the eigenvalues of Hi, see Eq. (JJ]). But it is easy 
to see that the large-scale diagonalization eigenvalues do not have this structure - we only 
obtain this in the limit iV max — > oo. Therefore we choose the eigenvectors corresponding to 
the N eigenvalues E n>m , n + m < iV max , thereby achieving this structure in the eigenvalues 
of H> s . 

In general, we wish to incorporate the symmetries of H' into the effective Hamiltonian 
H' cS . In this case, it was the separability and even eigenvalue spacing we wished to reproduce. 
In Sec. |V] we treat the two-dimensional Coulomb problem similarly. 

B. Numerical experiments with effective interactions 

The eigenvectors of the Hamiltonian H' differ from those of the the effective Hamiltonian 
H' efl . In this section, we first make a qualitative comparison between the ground states of 
each Hamiltonian. We then turn to a numerical study of the error in the CI method when 
using the effective interaction in a the model problem. 

Recall that the ground state eigenvectors are on the form 

oo 
n=0 

For H' eS , c n = for all n > iV max , so that the excluded space-part of the error concides 
with the excluded space-part of the exact ground state. In Fig. the coefficients c n for both 
H' and H' cS are displayed. The pointwise error is also plotted, and the largest values are 
seen to be around x — 0. This is expected since U(\/2x; A, 5) and the exact ground state is 
non-smooth there. Notice the slow spatial decay of the error, intuitively explained by the 
slow decay of the Coulomb interaction. 
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FIG. 5: Left: Plot of ground state coefficients of H' and H c g. Right: Pointwise error (in relative 
coordinate x) of effective Hamiltonian ground state ip e s(x) 

We now turn to a simulation of the full two-particle Hamiltonian (jSJ), and compare the 
decay of the ground state energy error with and without the effective interaction. Thus, we 
perform two simulations with Hamiltonians 

H = H' +v(xi) + v(x 2 ) 

= h(xt) + h(x 2 ) + v( Xl ) + v{x 2 ) + TVf ] 

and 

H cS = H' e{i + v{x x ) +v(x 2 ) 

= h(xi) + h(x 2 ) + v(xi) + v(x 2 ) + fV e gf\ 

respectively, where T is the centre-of-mass transformation, cf. Eq. ([7]). 

We remark that the new Hamiltonian matrix has the same structure as the original 
matrix. It is only the values of the interaction matrix elements that are changed. Hence, 
the new scheme has the same complexity as the CI method if we disregard the computation 
of Veg, which is a one-time calculation of low complexity. 

The results are striking: In Fig. Owe see that the ground state error decays as 0(A^"~^ X 57 ), 
compared to 0(iV~° x 95 ) for the original CI method. For iV max = 40, the CI relative error 
is AE/Eq w 2.6 • 10 -3 , while for the effective interaction approach AE/Eq fa 1.0 • 10~ 5 , a 
considerable gain. 
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FIG. 6: Ground state energy relative error for a two-particle simulation using the confinement 
potential V(x) = x 2 /2 + 4exp(— 2{x — 0.75) 2 ). For the CI method without effective interactions, 
we obtain a ~ —1.02, while the effective interactions gives a ~ —2.57. The electron density is 
superimposed on the potential plot. 

The ground state energy E used for computing the errors were computed using extrap- 
olation of the results. 

We comment that iV max ~ 40 is the practical limit on a single desktop computer for a 
two-dimensional two-particle simulation. Adding more particles further restricts this limit, 
emphasizing the importance of the gain achieved in the relative error. 

In a more systematical treatment, we computed the error decay coefficient a for a range 
of trap potentials x 2 + Aexp(— 2(x — /i) 2 ), where we vary A and /j to create single and double- 
well potentials. In most cases we could estimate a successfully. For low values of /i, i.e., 
near-symmetric wells, the parameter estimation was difficult in the effective interaction case 
due to very quick convergence of the energy. The CI calculations also converged quicker in 
this case. Intuitively this is so because the two electrons are far apart in this configuration. 

The results indicate that at iV max = 60 we have 

a = -0.96 ± 0.04 for H 

and 

a = -2.6 ± 0.2 for H eS 

for the chosen model. Here, 0.6 < /x < 1.8 and 2.9 < A < 4.7 and all the fits were successful. 
In Fig. [7] contour plots of the obtained results are shown. For the shown range, results were 
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FIG. 7: Estimates of a for CI calculations with (right) and without (left) effective interactions, 
unambiguous. 

These numerical results clearly indicate that the effective interaction approach will gain 
valuable numerical precision over the original CI method in general; in fact we have gained 
nearly two orders of magnitude in the decay rate of the error. 



V. DISCUSSION AND OUTLOOK 



A. Generalizations 



One-dimensional quantum dot models are of limited value in themselves. However, as 
claimed in the Introduction, the analysis and experiments performed in this article are valid 
also in higher-dimensional systems. 

Consider two particles in two dimensions. Let h(r) be the two-dimensional harmonic 
oscillator Hamiltonian (we omit the superscript in Eq. (j3J) for brevity), and let the quantum 
dot Hamiltonian be given by 

H = H' + v (n) +v(f 2 ), 

where 

H> = h{f x ) + h(r 2 ) + — X — . 

Wn - r 2 \\ 

The normalized centre-of-mass and relative coordinates are defined by 

- f i + f 2 _ r i - r 2 
K = ^— and r = ^— , 



23 



respectively, which gives 

H' = h(R) + h(r) + A 



y/2\\r\\' 

The h.o. eigenfunctions in polar coordinates are given by^ 

® n , m (r,e) oc e^ e rH L H (r 2 )e -, 2 /2 ; 

and the corresponding eigenvalues are 2n + \m\ + 1. Now, fT is further separable in polar 
coordinates, yielding a single radial eigenvalue equation to solve, analoguous to the single 
one-dimensional eigenvalue equation of Hi in Eq. (j3j). 
The eigenvalues of H' have the structure 



E n >,m',n,m = 2n' + \m'\ + 1 



where (pf, mf) and (n, to) are the centre-of-mass and relative coordinate quantum numbers, 
respectively. Again, the degeneracy structure and even spacing of the eigenvalues are de- 
stroyed in the CI approach, and we wish to regain it with the effective interaction. We then 
choose the eigenvectors corresponding to the quantum numbers 

2n' + \m'\ + 2n + m < N max 

to build our effective Hamiltonian H' cS . 

Let us also mention, that the exact eigenvectors ^ n ',m',n,m are non-smooth due to the 1/r- 
singularity of the Coulomb interaction. The approximation properties of the Hermite func- 
tions are then directly applicable as before, when we expand the eigenfunctions in h.o. basis 
functions. Hence, the configuration-interaction method will converge slowly also in the two- 
dimensional case. It is good reason to believe that effective interaction experiments will 
yield similarly positive results with respect to convergence improvement. 

Clearly, the above procedure is applicable to three-dimensional problems as well. The 

operator H' is separable and we obtain a single non-trivial radial equation, and thus we may 

apply our effective Hamiltonian procedure. The exact eigenvalues will have the structure 

3 

E nl>v>ml>n>l>m = 2n' + l' + - + e njl>m , 

on which we base the choice of the effective Hamiltonian eigenvectors as before. 

The effective interaction approach to the configuration-interaction calculations is easily 
extended to a many-particle problem, whose Hamiltonian is given by Eq. ([T]). The form of 
the Hamiltonian contains only interactions between pairs of particles, and V e s as defined in 
Sec. [TV] can simply replace these terms. 
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B. Outlook 



A theoretical understanding of the behavior of many-body systems is a great challenge 
and provides fundamental insights into quantum mechanical studies, as well as offering po- 
tential areas of applications. However, apart from some few analytically solvable problems, 
the typical absence of an exactly solvable contribution to the many-particle Hamiltonian 
means that we need reliable numerical many-body methods. These methods should allow 
for controlled expansions and provide a calculational scheme which accounts for succes- 
sive many-body corrections in a systematic way. Typical examples of popular many-body 
methods are coupled-cluster methods,- 1 ^ 1 ^ various types of Monte Carlo methods,— 
perturbative expansions,— 1 ^ Green's function methods,— ^2 the density-matrix renormaliza- 
tion group^ 1 ^ and large-scale diagonalization methods such as the CI method considered 
here. 

In a forthcoming article, we will apply the similarity transformed effective interaction 
theory to a two-dimensional system, and also extend the results to many-body situations. 
Application of other methods, such as coupled-cluster calculations, are also an interesting 
approach, and can give further refinements on the convergence, as well as gaining insight 
into the behaviour of the numerical methods in general. 

The study of this effective Hamiltonian is interesting from a many-body point of view: 
The effective two-body force is built from a two-particle system. The effective two-body 
interaction derived from an iV-body system, however, is not necessarly the same. Intuitively, 
one can think of the former approach as neglecting interactions and scattering between three 
or more two particles at a time. In nuclear physics, such three-body correlations are non- 
negligible and improve the convergence in terms of the number of harmonic oscillator shells.— 
Our hope is that such interactions are much less important for Coulomb systems. 

Moreover, as mentioned in the Introduction, accurate determination of eigenvalues is 
essential for simulations of quantum dots in the time domain. Armed with the accuracy 
provided by the effective interactions, we may commence interesting studies of quantum 
dots interacting with their environment. 
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C. Conclusion 

We have mathematically and numerically investigated the properties of the configuration- 
interaction method, or "exact diagonalization method" , by using results from the theory of 
Hermite series. The importance of the properties of the trap and interaction potentials is 
stressed: Non-smooth potentials severely hampers the numerical properties of the method, 
while smooth potentials yields exact results with reasonable computing resources. On the 
other hand, the h.o. basis is very well suited due to the symmetries under orthogonal coor- 
dinate changes. 

In our numerical experiments, we have demonstrated that for a simple one-dimensional 
quantum dot with a smooth trap, the use of similarity transformed effective interactions 
can significantly reduce the error in the configuration-interaction calculations due to the 
non-smooth interaction, while not increasing the complexity of the algorithm. This error 
reduction can be crucial for many-body simulations, for which the number of harmonic 
oscillator shells is very modest. 
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